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A continuous-time formulation of the Diffusion Monte Carlo method for lattice models is pre- 
sented. In its simplest version, without the explicit use of trial wavefunctions for importance sam- 
pling, the method is an excellent tool for investigating quantum lattice models in parameter regions 
close to generalized Rokhsar-Kivelson points. This is illustrated by showing results for the quantum 
dimer model on both triangular and square lattices. The potential energy of two test monomers as a 
function of their separation is computed at zero temperature. The existence of deconfined monomers 
in the triangular lattice is confirmed. The method allows also the study of dynamic monomers. A 
finite fraction of dynamic monomers is found to destroy the confined phase on the square lattice 
when the hopping parameter increases beyond a finite critical value. The phase boundary between 
the monomer confined and deconfined phases is obtained. 
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In the past decade substantial progress have been made 
in applying Monte Carlo simulations to study quantum 
lattice models. Yet there are still quantum models which 
are hard to simulate efficiently, even though they do not 
suffer from the infamous sign problem. Recently there 
has been a revived interest in the quantum dimer model 
(QDM) , a quantum model for which it is hard to engineer 
efficient Monte Carlo updates. 

The QDM was first proposed in the context of 
resonant valence bond (RVB) theories of high-Tc 
superconductivity0|. In the RVB theory^] pairs of spins 
form singlets which are approximated by dimers in the 
QDM. The RVB scenario is interesting as it implies a 
gapped spin liquid phase where holes can move freely, and 
provides an exotic mechanism for superconductivity 
Thus a particularly interesting feature of the QDM is 
the existence of a non-trivial point in parameter space, 
the Rokhsar-Kivelson (RK) point, where the model is 
exactly solvable ,1], and where monomers (holes) are de- 
confined. However it has been shown by exact studies on 
small lattices 1^ m, and repeated here for bigger system 
sizes, that the RK point is special and that it is sur- 
rounded by crystalline phases such that monomers are 
confined for generic values of the parameters in the QDM 
on the square lattice. This notion of a deconfined critical 
point has recently received attentionjS, iSj, and effective 
field theories describing the nature of small perturbations 
away from this point has been proposed. The situation 
on the triangular lattice is different, there it was shown 
recently^ that the spin liquid behavior exists also away 
from the RK point. 

While it is interesting to know how the QDM behaves 
when the dimers fully cover the lattice, it is even more 
important, at least with regards to superconductivity, to 
know how the model behaves when a finite fraction of 
monomers is introducedQ- The Monte Carlo method in- 
troduced here allows the study of this. We find that the 
crystalline monomer confining phase is destroyed when 



the product of the monomer concentration and hopping 
increases beyond a finite critical value. Our main re- 
sult besides the Monte Carlo algorithm itself is the phase 
boundary between the monomer confined and deconfined 
phases shown in Fig. 01 

The Hamiltonian of the QDM is 

H^-jJ2 (|ll)(=l + H.C.) +VJ2 (ill) ("I + !=)(= 

(1) 

where the summations are taken over all elementary pla- 
quettes of the lattice. The reason why the QDM is a 
difficult model to simulate stems from the fact that effi- 
cient updates of the classical dimer model are necessarily 
non-local. In an imaginary-time formalism a space-time 
configuration consists of sheets of dimer configurations 
where nearby sheets in the imaginary time direction may 
differ by the orientation of two parallel dimers at the 
same plaquette. While efficient non-local update tech- 
niques such as the directed-loop method^] can be uti- 
lized for the classical dimer models^^l' J^st updating the 
dimers within one space-time sheet, the QDM seems a 
formidable challenge in comparison, as all parallel sheets 
would need to be updated simultaneously. Because of 
this it is reasonable to expect better performance from 
Monte Carlo techniques that avoids the constraint im- 
posed by the periodicity in the imaginary-time direction. 

Diffusion Monte Carlo (DMC)jjl2,] is a genuine zero- 
temperature technique and contains no references to the 
periodicity in the imaginary-time direction. It is part of a 
general class of methods known as Projector Monte Carlo 
techniques, where a stochastic process is used to model 
the result of repeated matrix multiplications. DMC is the 
special case where the iterated matrix is exp(— iJAr). 
Thus it describes, for infinitesimal Ar, the imaginary- 
time evolution of the wavefunction. 

As shown by Henley^3| the dynamics of the QDM at 
the RK point can be obtained using classical continuous- 
time Monte Carlo. Thus a continuous-time formulation 
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of the DMC is needed. It is often stated in the ht- 
erature that DMC cannot be formulated in continuous 
(imaginary) time, and that repeated runs with decreas- 
ingly smaller time-intervals must be performed in order 
to quantify the error induced by a finite time-step. How- 
ever, as we will show, a continuous-time formulation of 
DMC is feasible. Moreover this method becomes identi- 
cal to a classical Monte Carlo at the RK point. The algo- 
rithm presented here do not work for spatially continuous 
systems and one should consult Ref. to reduce the 
time-step errors. 

The algorithm will be illustrated for a two-state sys- 
tem with Hamiltonian matrix elements: Hij = {i\H\j) — 
EnSij, where i,j = {1,2}, and Er is a reference en- 
ergy. The generalization to larger systems is trivial and 
is stated in Eq. I^. As in the conventional DMC one 
introduces M replicas (or walkers) each one representing 
a basis state of the system. The collection of replicas 
represent an instance of the ground state wavefunction 




where Mj is the number of replicas in state |j) (Mi + 
M2 = M). 

We will consider the time evolution operator for an 
infinitesimal time step dr. The action of the evolution 
operator on an instance of the state is 

/ M( \ / 1 - Hndr -Hi2dT \ ( 

\M^J \ -H2idT 1 - H22dT J \M2 )■ ^ ' 

We will now formulate a stochastic process that on av- 
erage yields the evolution equation above: In the time 
interval dr a replica in state \i) can undergo one out of 
four different processes: (1) with probability Pt(*) it can 
change state to j 7^ i, (2) with probability Poii) it 
can "die", that is Mi Mi — 1, (3) with probability 
PR{i) it can "replicate", that is an extra replica in state 
\i) is created, and finally (4) with probability Ps(i) it can 
stay unchanged in state As these are all possibilities, 

Prii) + Poi^) + Piii^) + Psi^) - 1, (4) 

and must hold for all states i = 1,2. Because the off- 
diagonal matrix element Hji is the only one responsible 
for transition between state i and j it is clear that 

Pt{i) = -H,,dT (5) 

where j ^ i. As usual on order to avoid the sign problem 
off-diagonal matrix elements are restricted to be negative, 
which will be assumed in the following. 

The increase in number of replicas in state \i) from 
processes acting on replicas in state |i) is 

m; - = [Pnii) - Prii) - Poii)] M,. (6) 



This implies when comparing to the diagonal elements of 
Eq. and using Eq. (0) that 

Poii) - PB.{i) = [Hu + Hj,) dT (7) 

where j ^ i. The right hand side of the above takes 
either a positive or a negative value. We choose Pr — 
whenever this value is positive and Pd = when it 
is negative. This choice implies that Pd and Pr are 
of the order dr as also holds for Pt- The probability 
conservation equation Eq. Q then implies 

Psii) - 1 - + H,, I - H,,) dT (8) 

Thus for most time intervals nothing happens to each 
replica. This is analogous to stochastic modeling of the 
radioactive decay problem, although here with three dif- 
ferent decay channels|lfil|. Thus we can simulate the 
imaginary-time evolution of one replica by generating ex- 
ponentially distributed decay times with decay constant 
\Hii ~\- Hji \ — Hji. Having obtained the decay time, the 
type of decay is determined stochastically proportionally 
to the respective probabilities Pt, Pd and Pr. 

The generalization to larger systems is immediate. In 
Eq. © Pr(«) should be changed to account for transi- 
tions to all possible states different from It follows 
that the only change to Eqs. Q and © is to set 

H,^^Y.HJ^■ (9) 

In an actual simulation each replica contains informa- 
tion about the state of the system as well as a "clock" 
indicating the starting time for the next evolution. The 
replicas are ordered in a list. They all start in the same 
state and with their clock set to 0. Each replica in the 
list is subsequently evolved up to a control time Tc, or 
until the replica dies in which case it is removed from 
the list. An evolution of a replica begins by generating a 
decay-time according to the exponential distribution. 
If Td > Tc the clock is set to Tc and evolution of the next 
replica starts. If however Td < Tc, the clock is set to Td 
and a random number is drawn to select the decay type. 
If the decay type is (1) the state of the rephca changes, 
if it is (2) the replica is removed from the list and if it is 
(3) a copy of the replica with clock set to Td is inserted at 
the end of the list. As long as the replica is not dead the 
evolution continues by picking a new decay time until Tc 
is reached. When the last replica in the list has evolved 
up to Tc, all replicas have the same clock-time, and mea- 
surements can be performed. The process is repeated by 
increasing Tc and starting over from the beginning of the 
replica list. 

The control times are included in order to perform pop- 
ulation control to avoid an explosion/implosion in the 
number of replicas. Population control is achieved by 
changing the value of the reference energy Er so as to 
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maintain a roughly constant number of replicas. While 
a constant Efi is innocuous, a time- varying one is not. 
Thus the simulation with population control is not sim- 
ulating exactly the time evolution operator. A way to 
correct this is to reweight the simulations by a factor 
that corrects for the time- varying reference energy {Ij. 

It is known that importance sampling reduces statisti- 
cal errors in DMC. Importance sampling is achieved by 
sampling the product of the wavefunction times a trial 
wavefunction instead of the wavefunction alone. The trial 
wavefunction is chosen such that branching due to the 
replicating and death processes is minimized. We note 
that according to Eq. branching is absent when the 
potential energy of a state equals the kinetic energy as- 
sociated with motion away from the same state. This is 
exactly the condition defining generalized RK points [l^. 
Thus no explicit trial wavefunction is needed at RK 
points. Formally the (implicit) trial wavefunction is sim- 
ply 1, the equal superposition of all basis states, which 
in fact is the exact ground-state wavefunction. Thus wc 
expect small statistical errors in the vicinity of RK points 
even without introducing explicit trial wavefunctions. 

We now turn to simulations of the QDM. First we cal- 
culate the columnar order parameter, Xcoi, defined by 

- i^((E-H(^-)(-ir^)' + (E-v«(-ir")') 

where nn (ny) is the number of horizontal (vertical) 
dimers belonging to the plaquette at f. r is an integer- 
valued coordinate labelling the centre of each plaquette 
on a lattice of size N = L x L. The sums are to be 
taken over all plaquettes. The results are shown in Fig.Q] 
and are carried out for significantly larger lattice sizes 
than in the exact studies Refs. In all simulations 

presented here M = 1000 replicas were used. Diago- 
nal observables were measured using the forward- walking 
techniaue|l6l| keeping histories for times typically of the 
order 100/J. Energies were measured using the standard 
growth estimator [1^. To correct for the time- varying Efj 
we kept track of Eji for typical times 200/ J and used the 
approach described in Ref. to reweight the simula- 
tion. 

Extrapolating the columnar order parameter to infinite 
lattice size we find a linear behavior in 1/L at the RK 
point extrapolating to zero within error bars, and 1/L^ 
corrections for V < J extrapolating to finite values. Our 
error bars and system sizes investigated strictly allows 
the conclusion that the order parameter is finite for V < 
0.99 J, although we believe it is finite for all V < J. 

The existence of a finite columnar order parameter 
implies that monomers are confined. In order to show 
this explicitly we also compute the energy of two static 
monomers at different separations. Two monomers with 
a horizontal separation x is inserted in a configuration 
with horizontal columnar order. The insertion causes the 
dimers directly between the monomers to be displaced 
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FIG. 1: Columnar order parameter as a function of V/J on 
square lattices for different linear lattice sizes //(labels on each 
curve). The inset shows the energy per site as a function of 
monomer separation for V/J = 0.9. Error bars are smaller 
than the symbol sizes. 
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FIG. 2: Energy per site as functions of monomer separation 
for different system sizes L(labels on each curve) on the tri- 
angular lattice. 

one lattice spacing. At V = J the energy is indepen- 
dent of separation, but for I^ < J the energy increases 
as the separation between monomers is increased. This 
is shown in the inset of Fig. for V/ J — 0.9 and implies 
the confinement of monomers. Naively one would expect 
the graphs to be symmetric around L/2 because of the 
periodic boundary condition. However for small square 
lattices with L even the conservation of winding numbers 
prohibits this symmetry. 

Fig. 13 shows the energy of two static monomers on 
the triangular lattice as a function of their separation for 
V/J^ 0.9. The contrast with the square lattice case is 
evident. On the triangular lattice the energy increases for 
a few lattice spacings and then drops slightly reaching a 
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FIG. 3: Columnar order parameter on square lattices with 
dynamic monomer concentration n = 1/32 and VjJ = 0.9 
as a function of 1/L for different values of t/J(labels on each 
line) . The lines are least square linear fits to the data points. 
The inset on the lower right shows energy versus separation of 
two static monomers in a background of dynamic monomers 
and dimers for n = 1/8, L = 16, V/J = 0.9 and tjj = 0.1. 



plateau. Thus there is no confining potential. This agrees 
with the conclusions reached in Ref. 9] . 

We now turn to the study of monomers with their own 
dynamics. Specifically the termjll 



dyn 
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H.c 



(11) 



is added to the Hamiltonian Eq. (Q. The sum is taken 
over triplet of sites i,j and k such that both i and j , and 
j and k are nearest neighbors. 

In Fig. O we plot the columnar order parameter as 
a function of 1/L for different values of t for a small 
monomer fraction away from the RK-point. For t > tc ~ 
0.2J the columnar order extrapolates linearly to zero as 
L ^ CO. For smaller values of t this order parameter 
is finite. Thus there exist a finite critical tc where the 
columnar order is destroyed. Increasing n to 1/16, and 
1/8 we find that tc drops inversely proportional to n. A 
computation of the energy of two static monomers in a 
background of dimers and dynamic monomers for t > tc 
is shown in the lower right inset of Fig. |21 This indi- 
cates that the dynamic monomers screen the interaction 
between monomers, rendering the potential deconfining. 
Calculating the values of tc for other values of V/J we 
obtain the phase boundary between the confined and de- 
confined phases shown in Fig. ^ 
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FIG. 4: The boundary between the monomer confined and 
deconfined phases in the doped quantum dimer model. 



The method presented here is applicable to any quan- 
tum lattice models that do not suffer from the sign 
problem. This includes quantum vertex models [l^ and 
bosonic models with internal degrees of freedom. 

Monte Carlo calculations were in part carried out using 
NorduGrid, a Nordic facility for Wide Area Computing 
and Data Handling. 
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